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7 3314831235296416541804 3350226 
79344499 14650587 22779760834101 
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This is Factorial 11000 


.» « 48453312 388261 330008827929 

056400152870915856 31099566 3618 

86 384 37 3297981078 37608 35911330 

025575841231852468 337956093952 
and 2748 zeros, 
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4 Worlds record 


...feally high precision calculation 


It all started with a tentative solution to Problem 
120 (from our issue number 36) which appeared in issue 85. 
The problem concerned the sequence of non-zero low-order 
digits of successive factorials. The first 1000 of those 
digits had been calculated, primarily as an exercise in the 
use of floating point BASIC. 


There then followed much discussion of this sequence of 
digits with David Ferguson, who established that the sequence 
does not cycle, but that repeating strings of any length can 
be found. 


We have reproduced the listing of the first 1000 
of these digits (taken from page 3 of issue 85), together 
with the next 2000 digits. Mr. Ferguson noticed that 
in the first 1000 digits, the sequence starting at the 
627th digit repeated the sequence starting at the 2nd digit, 
so that a repetitive cycle of 625 digits was evident even 
from the limited data that was given. 


(The tables are on pages 4 and 5) 


Not to doubt for a moment as emminent a worker as 
Ferguson, but in the interests of acquiring more data, it 
seemed expedient to extend the calculations. The program 
in BASIC would be much too slow; moreover, it could not be 
extended very far due to space limitations. How about 
rewriting the whole thing in 6502 machine language? 


This is quite feasible. To go to, say, factorial 
5000, it would be necessary to operate on at least 1300- 
Gigit numbers at any one stage. 
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The reasoning to arrive at that figure of 1300 digits 
is quite simple. To get to factorial 5000, concerning 
ourselves only with the low-order non-zero digits, it is 
necessary to anticipate how many low-order zeros there can 
be. In developing the factorials systematically, every 
time the argument is a multiple of 5, the function gains a 
low-order zero. Each time the argument is a multiple of 
25, the function gains two low-order zeros, and so on. Thus, 
for factorial 5000, we have: 


Multiples of 5: 1000 
Multiples of 25: 200 
Multiples of 125: 40 
Multiples of 625: 8 
Multiples of 3125: ee 


1249 


So if 1300 digits are carried throughout the calculations, 
and the low-order zeros are shifted off at the right hand 
end of the number, the 51 extra digits should provide enough 
protection to insure that the low-order, non-zero digit is 
correct. 


With the usual (lazy) technique of carrying one 
decimal digit per word, there is no problem of storage 
space; 1300 consecutive words of storage are cheap. 


Multiplication is aiways a problem on a machine like 
the 6502 which lacks a MULTIPLY operation, and even more so 
when the multiplier can get to 5000 and the word size of 8 
bits limits the natural range of numbers to plus or minus 
127. So the following scheme was adopted. 


The index of the factorial being developed, N, will be 
contained in 3 words of storage in a base-40 scheme, so that 
the equivalencies shown in Table T hold. 


For each new value of N tee 1234), there will be 


formed immediately in storage (again using 3 words for each 
number) the 9 possible multiples of N. Thus, to form 

the next value of (N!), it is necessary only to look up 

the proper multiple in this table for each digit of (N-1)! 


Text continues on page 6 ) 


24284484666264662648868288682662644484644846224 284 
484688682224 28224 2866264886826626444 8464484 6224 288 


8682662644484644846224284484688682224 2822428662646 
62642242888682886824484688682662644484644846224286 
62642242888682886824484666264 224 288868288682448462 
24 284484666264662648868288682662644484 644846224284 


484688682224 28224 286626466682662644484644846224288 
86826626444846448462242866264224288868288682448464 
48468868222428224 2866264224 28448466626466264886824 
484688682224 282242866264448468868222428224 28662646 
62642242888682886824484644 84688682224 2822428662642 
2428448466626466264886824484688682224 28224 28662644 
484688682224 28224 28662644484688682224 2822428662646 
6264.2242888682886824484688682662644484644846224286 
62642242888682886824484666264224 288868288682448468 


8682662644484644846224 28224 2844 8466626466264886826 
6264.22428886828868244846224 28448466626466264886822 
24284484666264662648868266264224 288868 2886824U846K 
48468868222428224 2866264224 28448466626466264886824 
48468868222428224 2866264448468868222428224 28662648 
868266 2644484644846224 28224284 48466626466264886826 
626422428886828868244846224 2844 2466626466264886822 
24284484666264662648868288682662644484 644846224282 


the low-order, non-zero digits of successive factorials. 


For the first 1000 factorials, the table at the upper left 


on the facing page is reproduced from issue 85. The table 
on this page contains the digits from N = 1001 through 

N = 2150, 50 digits per line. Triple digits are under- 
lined, to aid in tracing patterns. The table at the lower 
right on the facing page contains the digits for N = 2151 
through N = 3250. 


12642242888682886824484644846886822242822428662642 
2284484 66626466264886824 484 68868g224 2822428662644 
UB46BB682224 28224 28662648868266264484644B46224282 
24284484 666264 662648868266264 224 2BBB68288682448462 
nu 20uusugegesusee6use6d2zzi 284484 G6G26u66264886828 
3682662644484 64484 622h 28224284484 66626466264886826 
6264224 2888682886824484622u 28448466626466264886822 
24284N84666264662648968222428484665626466264886828 
868266264448464484622428448468868222428224 28662648 


8682662644484 644846224 288868266 2644484644846224284 
484688682224 28224 286626466264 224 288868288682448468 


B6B2662644N84644846224 2866264 224 2B886828B8682UU8465 
6264.2242888682886824484666264224 2888682B8682448464 
48468868222428224 2866264224 28448466626466264886824 
48468868222428224 28662644484 688682224 2822428662648 
8682662644484 64484 622428224 28448466626466264886826 
6264224 2898682886824484 622 28448466626466264886822 
24 28448466626466264886828868 26626444 84644846224282 
24284%84666264662648868266264224 288868288682448462 
242844846662646626488682224 284484666526466264886822 


The single, low-order, non-zero digits of the first 
thougand factorials. 


24284.484666264662648868266264224288868288682448462 
24284484666264662648868222428448466626466264886822 
2428448466626466264886828868266264448464U846224284 
48468868222428224286626488682662644484644846224288 
86826626444846448462242844846886822242822h 28662646 
6264.224288868288682448468868266 2644484644846224286 
62642242888682886824484666264224 28886828868 2448464 
484688682224 28224 286626466264224288868288682448468 
86826626444846448462242866264224288868288682448466 
62642242888682886824484622428448466626466264886828 
8682662644481 6448462242844846886822242822428662648 
8682662644484644846224 288868266264448U64U846224282 
24284484666264.66264886828868266264U48464484 6224284 


484688682224 28224 286626488682662644484 644846224288 
8682662644484644846224 288868266 2644484644846224282 
24284484666264662648868266 264224 288868288682448462 
242844846662646626488682224284484666264 66264886826 
6264224288868288682448464484 688682224 2822428662642 
2428448466626466264886824484688682224 2822428662644 
48468868222428224 2866264224 28448466626466264886828 
86826626444846448462242844846886822242822428662648 
8682662644484644846224288868266264448464U846224286 


N (decimal) N (base 40 scheme--expressed 
in hexadecimal notation) 


1 
2 
3 
Mm 
5 
6 
7 
8 
g 


For example, suppose we have the following sequence 
of digits in storage for factorial 1233: 


Nelo SSOC8... 


and we are in the process of developing factorial 1234. 
For each digit involved, there is a Carry to consider, and 
Carry is also 3 words of storage. 


The maximum value that Carry can have is almost nine 
times the maximum N, or 45000; in the notation we are using, 
Carry will appear as (03 O04 05) for the decimal number 4965. 
At the start of the cycle that develops a new factorial, 
Carry is set to (00 00 00). 


Then, working from right to left, to replace the 
digit 8, the algorithm calls for looking up the 8th multiple 
of N, which is (06 06 20), adding the Carry value at that 
point, and reducing the result to the same form. The 
digit corresponding to the unit's digit of that result is 
stored back where the 8 was, and what is left is the value 
of Carry to use for calculating the next digit to the left. 


Having produced a part of each of the factorials up 
to factorial 5000, it seemed feasible to go further and 
use the same techniques to develop complete factorials. 
The opportunity was at hand to break a world's record, 
which is always fun. 


One should not plunge into such a project without 
first considering how the result might be validated. The 
last significant work (reported in our issue number 3) was 
done by Timothy Croy in 1972, to calculate factorial 10000. 
Mr. Croy thoughtfully printed out many intermediate results, 
and these could be used as checks in the new calculation 
of factorial 11000. 


An independent calculation was made, to sum the 
common logarithms of the numbers from 2 to 11000; this sum 
(made on an SR-52) was: 


39680 .4999853 
which indicates that factorial 11000 has 39681 digits, 
the first few of which are 31621... Further, we can 


calculate, as we did before, the number of low-order zeros, 
which is 2748. 


A further check was furnished by Herman P. Robinson, 


using Stirling's formula, giving the first 45 digits of the 
expected result. 


There is one more interesting aspect to this problem. 
We can set up a work space of 39681 words, initialized to 
zeros except for a 2 in the low-order word, and start 
producing factorials with N = 3. Clearly, this will be 
extremely inefficient, inasmuch as the result for the first 
product is contained in one word, but the straightforward 
procedure will involve 39680 "multiplications" by zero. 
What is needed is some scheme to allow the multiplications 
to proceed only as far as is needed for each value of N. 
All sorts of schemes suggest themselves, but each one 
requires elaborate coding to implement it, and that code 
chews up execution time, too. Fortunately, there is a 
simple scheme that just happens to apply here. 


Since factorial 11000 has nearly 40000 digits, it will 
be safe to allow 4N digits for each factorial. Thus, the 
loop that forms each new factorial can be terminated at an 
address that is moved 4 to the left for each new value of N. 
This scheme is a little wasteful (for example, factorial 
1000 has 2568 digits, and this scheme is allowing 4000 
digits), but on the side of discretion. And, of course, 
altering a test address by 4 each time around the main loop 
requires very little code. 
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At that, the errors in the above scheme can be corrected 
if the whole procedure is halted every 1000 stages to record 
results and compare with previous values. 


For example, suppose a HALT is arranged at factorial 
1000, to check the leading digits with Croy's value (which, 
incidentally, was reproduced explicitly in our issue number 
3). With the machine halted, it is a simple matter to 
alter the test address (which is at this point in error by 
over 1400 words) and restart the procedure, heading now for 
factorial 2000, at which point another correction can be 
applied. 


When N is small, at the start of the calculation, 
this process runs at around four values of N per second on 
a 6502 microprocessor. At the other end, when N is 
approaching 11000, it runs around 36 seconds per N value. 
The total CPU time for factorial 11000 was 42 hours. Ite 
is estimated that nearly 39 billion instructions were 
executed. 


The final result is a world's record. Ig abs) ioe, 
as the saying goes, "a pivotal event in world history,” 
but it is a thoroughly satisfying job. Incidentally, in 
one sense, the result is the largest number ever calculated. 


There was no attempt made in writing the program for 
factorial 11000 to conserve storage space or CPU time. 
With only modest effort, it would be feasible to pack two 
decimal digits in each 6502 word, this saving storage space 
(indeed, the 6502 permits decimal arithmetic directly for 
such an arrangement). 


And consider those 2700 trailing zeros. From time 
to time during the calculation (say, every 1000 stages), 
all the trailing zeros could be shifted off. There must 
be better schemes for performing the multiplications, and 
more efficient ways to adjust each multiplication to its 
proper length. All the calculations can be done on any 
microprocessor. 


So: it should not be too long before we have 


factorial 12000. 


&) 
O 


O 
O 


Fred Gruenberger 


“OT Jo azemod agedoud 3ay4 pue 
S4TSTP queotsitusts oG 4Ssity |ayq BupMoys 


fgouenbas JFooeucgdTy 34uq JO sulla} paqoatas 


OTTOSLTISSESQHELETCLLILHETEOLEEESSOnHE9S6ETHShSOTO’ 2 
79917800R8LOgLrEegcTt 7lS6r9t4LyLOOBOTSEZeSZeE TLOnSESere' T 


OLLEBOS9OTREZTLOE LEHE6eLSOLESHEE LOSTORBOHERRT69S6eS2 "2 
nOLOZEISROHIEELLONRIhHeTLEESTELS9RSLEBSOTSBE OF EE * T 


SHE TE9ZZRERcE HEGOS6OTS HOSLE OH 6ETE HEE 69995656965 6Sh2S ° 2 
EOSETIORST HE HSETTLOSREE 19 hhHEBEZO0S HSSLT H692T98Z09S" T 


ZTE 6SESEOEST HSHEZSLOEE GE LOBERRESHSeHGb6ERg90T9nELI’ 2 
St90SCZ66LSSETE LET B9RLERE te prOR69L SOQEETHLLT LOTS’ T 


BEQILESSORSZSB0EEE B2CSSOOE SIRE 9STZESR6LEQ9E HS90T6’ 2 
69ETETLOOSHOSLHOTITLLLEQE TS6LT LOG6E9ZEHZE SHCEE BRRbL* T 


GOTE90QTTB0E SESH LEZHEOT OE 6R76E 2S H66TLSHegHTQOETTgO’€ 
EEZO9TH?PZOHHESTOSGLSLTILYB9LE00SSHLEE HeOHE LOGE HeHO6’ T 


TTLOL9NS6099E LOG LEQROTEEESESLETE0STLICE H6L18E 6225S" € 
T2OECSESTOS6HEEZ6OCLT Ent 6TELOGLZTOE OE LISSSE GE HS6T Zz 


G6eS 10669222 LLE9HSLEOE HIG TOLEHTLS99R8TZOeTLSeet7eS9’ E 
LUTSLOSectHE HLOOBHE 6QELLLOZIBHTLTSS9EEOLHGE LOTERSC 2 


O6tHTHE0E LOH LyLLETE9GELETOTOE H6HO60E 9HSGe2HESOTSEOGL’ 
CO6LTEGLOE HIGEZE TSE QQ 609HnT Bb6SeLlrHETgge LOT2rLHeOHZEe “2 


99LOLOEE 669066E TE B90029L69£ S8S6299ST9L2Z6E 9506898" £ 
ZEELONOS H9EOLH LON? HeTERQZE T6655 67S 6n08Q0£ HETSTTILOE *Z 


OBTSEESE 9R6SREE R1 Ez6S9rE hLEL60LE 6SQ0L] ngct On6G0R6’ € 
SSLLOLHLONTecERTZELOTETG0R90£E 61099 LLZTLLOVeHTO9n’ Zz 


. Fibonacci.. 


A much simpler high precision problem is that of 
calculating large terms of the Fibonacci sequence: 


See tere es 7 O83) 2), * 34,2. 


Tables of pairs of terms in this sequence have appeared in 
issues 25, 30, and 34, going as high as terms 20000 and 20001, 
at which point we were dealing with numbers of 4180 digits. 


The scheme outlined in Flowcharts U and V is simple 
and straightforward. As with the factorial calculation, no 
attempt was made to conserve either storage space or CPU time. 


The Fibonacci sequence gains one decimal digit about 
every 4.78 terms. Again, rather than devising an elaborate 
algorithm for controlling the number of digits being operated 
on, the simple rule was adopted of increasing the field 
length of each number by one digit every 4 terms. 


The accompanying table gives some new results, which are 
truncated from the exact explicit values; that is, they are 
NOT rounded from the 5lst digit (as they were in the table 
in issue 34). 


With access to 48K bytes of storage, even the one- 
digit-per-word scheme could extend this table quite far. 
For the effort of packing two digits per byte, a truly giant 
step could be taken in extending this table. 


The ratio of successive terms, F, of the Fibonacci 
sequence approaches: 


MSE 
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Herman P. Robinson writes “ 


-..the error in t obtained by 
dividing Fyo991 bY Fyqooo Will occur in the 16719th decimal. 
If that digit is greater than 6, the error could carry over 


to the 16718th decimal." 
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Givoli’s Problem 


On the facing page are two tables. In Table W, the 
middle column shows the consecutive prime numbers. The 
first column indexes them, taking 2 as the first prime. 

The column headed S is the progressive total of the primes. 


Nahman Givoli, of Tel Aviv, introduces this problem: 
If (p;) is the sequence of primes, 
2, 3, 5, 7, «++, then for each natural 
number d, find the smallest number N(d) 
such that d divides the sum of the 
first N primes (p,). 


Table Y shows the start of the desired table. Each 
of the natural numbers, d, is listed, together with the 
number, N, of the first prime for which d divides S. 


Givoli raises a question about this distribution: 
Can it be proved that an N(d) exists for every a? There 
is no reason to suppose otherwise, and yet for some small 
values of d the respective N values are very hard to find. 
This appears to be especially true for d's which are 
multiples of 6 


The mention of 6 in connection with primes rings a 
bell, since 6 appears to be the most frequent difference 
between successive prime numbers. However, Givoli's 
conjecture does not hold good consistently; we have these 
extensions of Table Y: 


60 103 25800 
66 175 83292 
70 209 122920 
72 119 35568 


and the last values of d for which we could not find an N 
were 303 and 326, neither of which is a multiple of 6. 
However, the problem warrants extensive investigation. 
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Three dice are tossed. The sum of the uppermost 
spots will be a number between 3 and 18. Out of 216 
tosses of three dice, the sums have the expected 
distribution as follows: 


Four tosses are made of the three dice, to locate 
one cell in each of the four arrays P, Q, R, and S. For 
the selected cells, the centers from P to S are connected 
by a straight line, and similarly for the selected cells 
in Q and R,. 


These two straight lines intersect within square ABCD 
and possibly within square EFGH. The intersection could 
lie exactly on the border of the inner square EFGH, and 
in that case is to be discarded. 


For example, if the dice call for the selection of 
cell 16 in array P, 18 in S, 3 in Q, and 11 in R, then 
the intersection of the two lines wili fall on the border 
of the inner square. 


The Problem, then, is: in the long run, what 
fraction of the intersections will fall inside EFGH? 


rs is [| 
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Exploring Random Behavior -- 5 


So you think you understand how random processes will behave? 


Problem Solution 


Dr. Rudi Borth, of Toronto, offers a different 
algorithm for the TAKE/SKIP7 problem that appeared in 
issue 88. The problem was this: 


Start with the positive integers. For the first 
level of sieving, Take the first two numbers, 

Skip the next 3, Take the next 4, Skip the next 5, 
Take the next 6, Skip the next 7, and so on, 
indefinitely. 


Each lower level operates in much the same way. 
Level K starts by Taking (K+1), Skipping (K+2), 
Taking (K+3), Skipping (K+4), and so on, 
indefinitely. 


What numbers will survive all levels of this sieve? 


What follows is from Dr. Borth's letter: 


After the first two terms (1,2), each level K = L 
adds one number to the series; namely, the number in position 
P=K+1, To determine its value, all we need is to 
trace the number back to its position in level L = 1 (the 
natural numbers) where position equals value. 


A recurrence relation which computes position Q of a 
number in level (L - 1)Sfrom its position P in level L is 
obtained from the fact that the difference Q@ - P equals the 
sum of the I "Skips": 

Moe Ga) + (its jet <6: & (LtlI-1) =*Ibe+ 1° 


which have produced level L from level L-1. Thus, 
Q=P + I(L4+I). 


The required I is determined from the fact that I is 
also the number of "Takes" which have created level L; that 
is, the greatest possible number of terms in the sum: 

Meee) FG) + oot (Lite a(l - 1)) = a¢Lti-1) 


such that the sum does not exceed position P - l. Jie ate 
is the positive root satisfying the second-degree equation 


Woe ce a0 Sali ie S A 
then the required I = INT(I'). 


An implementation of this algorithm on an H-P 9825A 
took 11 seconds to produce the first 30 members of the 
series. It seems that the time for each cycle is proportion- 
al to K rather than ok, 


X--X--K--X-- KK KK KK KH K--K--X 


Dr. Borth is, of course, quite correct, and his 
solution to the problem is delightful. The presentation 
that appeared in issue 88 can be justified on these 
grounds, however: 


(1) The overall objective was to show the great 
difference in execution time between floating point BASIC 
and machine language, using the identical algorithm in 
both situations. For this purpose, it is irrelevant how 
inefficient the solution is. 


(2) The suggested problem offered another opportunity 
to demonstrate the "bucket brigade" method. 


(3) We can't all devise brilliant and clever solutions; 
most of us just plug along. Thank heaven our computers are 
so fast. 


A flowchart for Borth's 
algorithm is given on page 20. 


John W. Wrench has written to point out that the result given 
in issue 89, page 20, is not correct beyond the 23rd significant 
digit. A better result (calculated by Herman P. Robinson) is: 


8.70003 66252 08194 50322 24098 59113 00497 11932 97949 74289 20921 ... 


Mr. Wrench also points out that proper terminology calls for 

labelling this result the limit of R, 48 n approaches infinity. 
Finally, the last result given on page 20 of issue 89 seems to be oo 
exactly half its correct size, and should be 8.69866568. od 
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